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Abstract 

In this short note we review deterministic simulation of biochemical 
pathways, i.e. networks of biochemical reactions obeying the law of 
mass action. It is meant as a basis for the MATLAB code, written 
by the author, which permits easy input and simulation of general 
biochemical networks. This work was carried out for the European 
Project 'CardioWorkBench'. 

1 Introduction 

These notes give a short guidance to deterministic simulation of biochemical 
pathways. We follow the approach introduced by Ullah et al. [5] , basis of the 
system biology sbtoolbox [4] freely available from [1]. For a comprehensive 
list of system biology software see [2]. 

The aim is to easily input and simulate a general biochemical network 
of reactions obeying the law of mass action. The pathway specifications 
(rates, reactants, and products of each reaction) have to be easily and flexibly 
importable. And the evaluation of the resulting system of ODEs expressing 
the model's dynamics have to be computationally efficient. We optimise 
previous implementations and show how to efficiently calculate the ODEs 
jacobian. We also mention how to easily include parameters like phases and 
volume scalings. A user-friendly MATLAB code implementing the algorithm 
discussed in these notes can be downloaded from [3]. 
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2 Modeling 



We consider simulating a biochemical reaction pathway (network) assuming 
that it can be decomposed into unidirectional elementary reactions and that 
the law of mass action can be applied to each reaction. 

Let X = (X Tt )ivxi represent the molecular species, R = (R m )i X M, the el- 
ementary reactions, k = (k m )i X M the rate coefficients. The set of elementary 
reactions forming the pathway can be written as 

N N 
Rrn ■ ^ ^ t>nm,X n ^ ^ ^ r nm X n , 171 1, . . . , M, 

n=l n=l 

where the reactions' stoichiometric coefficients l nm and r nm are non-negative 
integers. The step change in the number of molecules X n due to reaction R rn 
is given by 

dnrn Turn ^nrn- 

We collect these coefficients into the input and output stoichiometric ma- 
trices L = (l n m)NxM and R = (r nm ) NxM , and the step change matrix 

D = (^nm)jVxM- 

In practice, it is convenient to input the information contained in the 
matrices L and R by storing their (few) non-zero entries values and indices, 
as it is done in any matrix sparse representation. This process is described 
below with an example. 

The law of mass action implies that the dynamics of the molar concen- 
trations x n = [X n ], n = 1, . . . ,N is described by a system of ODEs, each 
called a rate equation, that we can easily write down in terms of the reac- 
tions' coefficients. The set of molar concentrations x = (x„)at x i must satisfy 
the system of first order ODEs 

x = f(x) with f = (/„)jvxi and (1) 

M / N \ 

/„(x) = ^6U fc m JI^ im • (2) 

m=l \ i=l / 

The parameter d nm is different from zero if the n-th variable is involved in the 
m-th reaction, the rate of which is given by the expression in parentheses. 
Note that all this can be generalized to include more complex dynamics 
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by allowing the stoichiometric coefficients to be non-negative real numbers 
(Generalized Mass Action), see [6]. 



Example. Consider the basic enzyme-kinetic reaction 



X\ + X2 X3 — >■ X\ + X4. 

k 2 

We decompose the reaction pathway into the elementary reactions 

fci. 



R\ : Xi + A 2 

k 2 



R? ■ X* 



R3 '■ A3 



A 3 

x 1 + x 2 

Xx + A 4 
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To this sequence of reactions we associate the following coefficient matrices 
used to compute (2): 



D 



These matrices can be obtained from a compact representation, that facili- 
tates the input of the model data. For instance, the matrix L specifying the 
reactants, can be computed from 



L 1 



The first row of L 1 (the superscript "i" stands for index) tells us that the 
first reaction involves the first and second reactant, and so on. The matrix 
L v ("v" for value) collects the related non-zero stoichiometric coefficients 
l nm . All remaining entries are filled with zeros. Similarly, we decompose the 
matrix R which contains the information about the reactions' products into 
two matrices R* and R v . 

In realistic biochemical pathways, there are many reactants and reac- 
tions involved but few reactants acting on the single elementary reaction. 
Consequently, most of the entries of the matrices L, R (and D) are zero. 
The compact (sparse) representation above represents a good compromise 
between user friendliness and computational efficiency. 
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3 Implementation and extensions 



The process of simulating a biochemical pathway can be decomposed into 
two tasks: the first is to set up the ODE system by obtaining the matrices 
of stoichiometric coefficients L and R, and the vector of reaction rates k; the 
second is the actual numerical solution (time-stepping) of the system with the 
given concentrations initial values. The computational cost of time-stepping 
is dominated by the evaluation of the ODEs' function, thus it is crucial that 
this is done efficiently. 

The implementation of f(x) presented in [5] uses the matrices L and D. 
Letting xx = (x, . . . , x)tv x m be a matrix whose M-columns are copies 1 of x, 
the expression in (2) can be evaluated as follows: 

f(x) = D* (k.*prod(xx. A L,l))', 

where * denotes matrix multiplication, .* and . A entry- wise matrix multi- 
plication and exponentiation, prod(-,l) multiplication along columns, and 
' transposition (these notational conventions are in accordance with MAT- 
LAB 's syntax). 

In [5], it is also proposed a code that speeds- up the evaluation of f(x) by 
avoiding multiplications by zero entries. This can be achieved more simply 
by using the compact representation above, and, if sparse matrices represen- 
tation is implemented, by creating D as a sparse matrix. A further straight- 
forward computational optimization is obtained by multiplying each column 
of D with the corresponding rate constant in k. Denoting the result of such 
multiplication by Dk, the computation of f(x) is reduced to the following 
steps: 



1. using L 1 , form the matrix X 



X\ £3 £3 

x 2 



2. evaluate f (x) as: f (x) = Dk * prod(X . A L V , 1)', 

with the convention that the empty product 0° = 1. 

Biochemical models often include parameters like, for instance, phases 
and volume scalings. It is useful to implement such parameters separately 
from the stoichiometric coefficients. This is easily done by multiplying the 
appropriate entry of Dk by the given parameter. 



1 in MATLAB, this is achieved by the command xx=repmat(x,l,M) 
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As the ODE system generated by biochemical pathways is non-linear, it 
is also useful to have a routine that evaluates the jacobian Jf (x) of the ODE 
system function f . The entries of the jacobian are given by: 

9fn (x) = f; dk nm (l jm x^- lfi) J] x l A , j,n = 1, . . . , N. (3) 



We proceed with the evaluation of the jacobian processing it by columns, 
i.e. evaluating for j — 1, . . . N. 

Let j G {1, . . . , N} be given. The evaluation of (3) is simplified by limiting 
the summation to the reactions involving Xj, and by limiting the multiplica- 
tions to the reactants involved in the given reaction. This latter simplification 
is already embedded in our compact representation. As for the summation, 
we define reduced matrices Dk^, X Xj , and obtained from the correspond- 
ing matrices Dk, X, and L v by considering only those columns associated to 
reactions that have Xj among the input reactants. 

We need to derive the entries of X Xj . A L v with respect to the variable Xj. 
To this end, we define a new row-vector l Xj collecting the values lj m with 
m such that reactant of reaction m, and a new matrix LZ obtained 

J Xj 

from L w x . by replacing the entries indexed jm with lj m — 1. Notice that all 
such matrix manipulations are made particularly easy by MATLAB built-in 
vector/matrix manipulation tools. In this way, we can express (and calculate) 
the j'-th column of the jacobian matrix in compact form as: 



df 

dxj 



x) Dk,, ^!, .* P nKl(X, . 14 ))'. 



Finally, let us discuss the problem of numerically solving (1). The main 
characteristics of such system of ODEs are the following. Unless all elemen- 
tary reactions are of 0-th or 1-st order, the system is nonlinear. Moreover, 
the system is generally stiff. Thus, it is compulsory to consider variable 
time-stepping and employ stiff ODE solvers. Standard ODE solvers pack- 
ages include robust stiff solvers (for instance, MATLAB's odel5s), that are 
fast enough if only a few simulations are needed. If, on the other hand, a large 
number of simulations is required as in parameter estimation, then it may 
be preferable to code ad hoc solvers that take into account the peculiarities 
of the ODE systems generated by biochemical pathways. 
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